function model_output=gft_ces(V_k,lambda,beta,gamma,theta,rho,K)
%This function solves for the gains from trade in the CES case

%Solve a single equation for reference sector lhat, then just use formulas
    function exdem1=exdem1_calc(V_k,lambda,beta,gamma,theta,rho,lhat1)
        term1=lambda.^(-(1-rho)./theta).*beta./V_k;
        term2=lambda(1).^(-(1-rho)./theta(1)).*beta(1)./V_k(1).*lhat1^(1+(1-rho)*gamma(1));
        term3=(term1/term2).^(1./(1+(rho-1).*gamma)).*V_k;
        exdem1=sum(term3)-1;
    end
%Initial values and options
lhat1_init=1;
options=optimset('TolX',.0025,'TolFun',.0025,'MaxFunEvals',500);
%Solving the model
fun=@(lhat1)exdem1_calc(V_k,lambda,beta,gamma,theta,rho,lhat1);
[x fval]=fsolve(fun,lhat1_init,options);
lhat1=x;

%Now plug into the formula for general lhat
t1=lambda.^(-(1-rho)./theta).*beta./V_k;
t2=lambda(1).^(-(1-rho)./theta(1)).*beta(1)./V_k(1).*lhat1^(1+(1-rho)*gamma(1));
lhat=(t1/t2).^(1./(1+(rho-1).*gamma));

%Now plug into utility formula
uhat = sum((lhat.^(gamma).*lambda.^(-1./theta)).^(1-rho).*beta)^(1/(rho-1));
model_output.welfare=1-uhat;

end